GravitationalWaves: a Python package designed to simulate, detect, and analyze continuous gravitational wave-forms. The package also takes in observed data for comparison in detection and analysis.
This report details through the software design decisions, usage
guidelines, and extensibility of the GravitationalWaves package system,
with respect to existing software libraries. The software libraries used
for the basis of this project include legwork and
riroriro. Legwork is a python package for computing
evolution and detection rates of gravitational-wave sources discovered
by space-based detectors such as LISA. Riroriro is a python package for
simulating and evaluating gravitational waves.
Gravitational waves, first theorized by Einstein, arise during events such as the collision of supermassive binary black hole mergers. A binary black hole is a system consisting of two black holes in close orbit around each other. These collisions are so powerful that they create distortions in space time, known as gravitational waves.
Definition. Gravitational waves are invisible distortions in space time, caused by massive events such as collisions between two black holes or neutron stars.
The Laser Interferometer Space Antenna is a proposed space probe to detect and accurately measure gravitational waves—tiny ripples in the fabric of spacetime—from astronomical sources. The space-based gravitational wave detector LISA presents a new view of gravitational waves by focusing on lower frequencies \(\left( 10^{−5} \lt f / \mathrm{Hz} \lt 10^{-1} \right)\) than ground-based detectors. This enables the study of many new source classes including mergers of supermassive black holes, extreme mass ratio inspirals, and cosmological GW backgrounds.
The importance of studying gravitational wave-forms stems from the idea of detecting and using gravity to estimate other dynamical astrophysical phenomena, giving enormous potential for discovering parts of the universe that are invisible to the eye, such as black holes and other unknowns.
GravitationalWaves is a Python package for simulating
the gravitational wave-forms of binary mergers of black holes and
neutron stars, computing several properties of these mergers and
wave-forms, and evaluating their detectability. In addition, the package
also takes in observed data from gravitational wave detectors, such as
the Laser Interferometer Space Antenna (LISA), to compare data and
predict detection rates for stellar-origin binaries. Specifically, we
implement the expressions for the evolution of binary orbits due to the
emission of gravitational waves, equations for the strain amplitudes and
signal-to-noise ratios of binaries, and approximations for the LISA
sensitivity curves.
Similar to the existing software libraries of legwork
and riroriro, the GravitationalWaves package represents an
implementation of gravitational wave tools, focusing on a mathematical
approach for analyzing different sources. Just like legwork, this
package provides mechanisms to calculate orbits strain amplitudes,
compute sensitivity curves, and calculate signal to noise ratios of a
collection of gravitational wave sources. Furthermore, the package
provides tools and plotting functions for visualizing such results. Just
like riroriro, this package provides gravitational waveform simulation
that’s based on a computational implementation of an earlier theoretical
waveform model that uses post-Newtonian expansions and simplified
Einstein field equations.
GravitationalWaves uses simulated data for gravitational waves and observations from detectors such as the Laser Interferometer Space Antenna (LISA).
Simulated Data: We obtain the simulated data by implementing the package functions for constructing simulations of gravitational waveform signals from a binary merger of two black holes or two neutron stars and outputting the data of such signals in terms of frequency and strain amplitude.
i_amp = ins.inspiral_strain_amplitude(i_Aorth, i_Adiag)
i_freq = mat.frequency_SI_units(i_m_omega, M)| Strain Amplitude | Frequency |
|---|---|
| 7.9097e-22 | 10.0066 |
| 7.9105e-22 | 10.0081 |
| 7.9112e-22 | 10.0096 |
| 7.9121e-22 | 10.0111 |
Observed Data: The observed data contains LISA-validated binaries obtained from the Kupfer et al. 2018 study. These observations are preloaded to the package and used to help generate the Source class based on these provided binaries.
| Source | f (mHz) | \(\beta\) (mas) | \(\sigma_\beta\) (mas) | d (pc) | \(\sigma_d\) (pc) | \(\mathcal{A}\) | SNR |
|---|---|---|---|---|---|---|---|
| AM CVn type | |||||||
| HM Cnc | 6.22 | NA | NA | [5000] | 6.4 | 211.1 ± 3.18 | |
| V407 Vul | 3.51 | 0.095 | 0.327 | 1786 | 667 | 11.0 ± 5.9 | 169.7 ± 2.17 |
| ES Cet | 3.22 | 0.596 | 0.108 | 1584 | 291 | 10.7 ± 4.6 | 154.3 ± 2.09 |
| SDSSJ135154.46-064309.0 | 2.12 | 0.596 | 0.313 | 1317 | 531 | 6.2 ± 3.5 | 21.8 ± 0.24 |
| AMCVn | 1.94 | 3.351 | 0.045 | 299 | 4 | 28.3 ± 3.2 | 101.2 ± 0.96 |
| SDSS J190817.07+394036.4 | 1.84 | 0.954 | 0.046 | 1044 | 51 | 6.1 ± 2.4 | 20.3 ± 0.13 |
| HP Lib | 1.81 | 3.622 | 0.052 | 276 | 4 | 17.5 ± 3.9 | 43.7 ± 0.28 |
| PTFl J191905.19+481506.2 | 1.48 | 0.550 | 0.327 | 1338 | 555 | 3.2 ± 1.8 | 4.0 ± 0.02 |
| CXOGBS Jl 75107.6-294037 | 1.45 | 1.016 | 0.146 | 971 | 156 | 4.2 ± 1.8 | 4.5 ± 0.02 |
| CR Boo | 1.36 | NA | NA | 337^a | ±44^a-35 | 13.4 ± 4.2 | 21.9 ± 0.13 |
| V803 Cen | 1.25 | NA | NA | 347^a | ±32^a-27 | 16.0 ± 5.4 | 26.2 ± 0.17 |
| Detached white dwarfs | |||||||
| SDSS J065133.34+284423.4 | 2.61 | 1.000 | 0.476 | 933 | 493 | 16.2 ± 8.6 | 90.1 ± 1.13 |
| SDSS J093506.92+441107.0 | 1.68 | NA | NA | 645^b | 41^b | 29.9 ± 7.7 | 44.9 ± 0.31 |
| SDSS J163030.58+423305.7 | 0.84 | 0.937 | 0.270 | 1019 | 357 | 11.5 ± 4.9 | 4.6 ± 0.03 |
| SDSS J092345.59+ 302805.0 | 0.51 | 3.340 | 0.173 | 299 | 10 | 26.4 ± 6.5 | 5.6 ± 0.06 |
| Hot subdwarf binaries | |||||||
| CD-30°11223 | 0.47 | 2.963 | 0.080 | 337 | 9 | 41.5 ± 1.8 | 4.9 ± 0.04 |
Limitations: Both simulated and observed calculations apply a low-order post-Newtonian description of gravitational wave emission, meaning that the data do not account for higher-order effects.
The following table presents the fundamental package dependencies:
| Package | Version | Type | Description |
|---|---|---|---|
| python | \(\geq\) 3.7.0 | Programming language | High-level, general-purpose programming language |
| pip | \(\geq\) 21.0.0 | Python package | Recommended package installer for Python |
| numba | \(\geq\) 0.5.0 | numba | NumPy-aware optimizing compiler using LLVM |
| numpy | \(\geq\) 1.16.0 | numpy | Fundamental package for array computing |
| astropy | \(\geq\) 4.0.0 | astropy | Astronomy and astrophysics core library |
| scipy | \(\geq\) 1.5.0 | scipy | Software for mathematics, science, and engineering |
| matplotlib | \(\geq\) 3.3.2 | matplotlib | Library for creating various visualizations in Python |
| seaborn | \(\geq\) 0.11.1 | seaborn | Library for making statistical graphics in Python |
| schwimm | \(\geq\) 0.3.2 | schwimm | Common interface for parallel processing pools |
| legwork | \(\geq\) 0.2.4 | legwork | Package for analyzing gravitational waves |
The structure of the GravitationalWaves github repository is as follows:
From the above, we can see the top-level package called
GravitationalWaves, with six modules:
strain.py, psd.py, utils.py,
visualization.py, source.py,
wavesim.py. There’s also the subpackage called
tests located in the top-level package. There’s also an
examples and docs directory that contains demos and documentation of the
package, though, they’re not packages since both don’t include an
__init__.py.
The software components for the GravitationalWaves
package include the package’s modules, classes, and functions. A Python
package is any directory with an __init__.py file, and this
file gathers all package-wide definitions to import every module in the
package. Python modules are an essential abstraction layer, which allow
separating code into parts holding related data and functionality. For
example, a layer of a project might handle interfacing with user
actions, while another handles data manipulation. To separate these two
layers, we’d regroup all user interfacing functionality into one file
and all data manipulation functionality into another file.
For the following component specification, we define each components’ metadata, interface, and implementation. The metadata consists of a component’s name and description. The interface specifies the list of inputs/outputs and their properties. And the implementation specifies how to execute the component instance, pass inputs to the component code, and get the component’s outputs. After the component specification, I introduce the component code, which implements the logic needed to perform a step in the package workflow.
Name: Source (class)
Description: The source module provides a simple interface to the functions in the remaining modules. This module contains the Source class for analyzing a generic set of gravitational wave sources, stationary/evolving and circular/eccentric.
Annotations:
"""
Class for generic gravitational wave (GW) sources.
"""
Inputs:
m_1,
description: Primary mass,
type: Int, optional:
False}m_2,
description: Secondary mass,
type: Int, optional:
False}ecc,
description: Initial eccentricity,
type: Int, optional:
False}dist,
description:
Luminosity distance to source, type:
Int, optional: False}n_proc,
description:
Number of processors to split, type:
Int, default: 1,
optional: True}f_orb,
description: Orbital frequency,
type: Int, optional:
False}a,
description: Semi-major axis,
type: Int, default:
None, optional: True}position,
description: Sky position of source,
type: Int, default:
None, optional: True}weights,
description:
Statistic weights for each sample, type:
Int, default: None,
optional: True}gw_lum_tol,
description:
Allowed error on GW luminosity, type:
Int, default: 0.05,
optional: True}stat_tol_e,
description:
Frequency change above which a binary is stationary,
type: Int, default:
1e-2, optional: True}interpolate_g,
description:
Whether to interpolate g(n,e), type:
Int, default: True,
optional: True}interpolate_sc,
description:
Whether to interpolate LISA sensitivity curve,
type: Int, default:
True, optional: True}Outputs:
snr, type:
Float, description:
Signal-to-noise ratio}Implementation:
import GravitationalWaves.source
import astropy.units as u
source1 = source.Source(m_1 = 11 * u.Msun, m_2 = 11 * u.Msun, ecc = 0.3, dist = 9 * u.kpc, f_orb = 1e-4 * u.Hz, interpolate_g = False)Name: strain
Description: A program that uses functions from
legwork to calculate signal-to-noise ratio of sources based
on four different cases. For the different cases, the module contains
four functions that cover whether a source is circular-stationary,
circular-evolving, eccentric-stationary, or eccentric-evolving.
Annotations:
"""
Computes strain amplitude at `n` harmonics
"""Inputs:
m_c,
description: Chirp mass of each binary,
type: Float, optional:
False}f_orb,
description:
Orbital frequency of each binary, type:
Float, optional: False}ecc,
description: Eccentricity of each binary,
type: Float, optional:
False}n,
description:
Harmonic(s) to calculate the strain, type:
Int, optional: False}dist,
description: Distance to each binary,
type: Float, optional:
False}position,
description: Sky position of source,
default: None, type:
Array, optional: True}interpolated_g,
description: G(n,e) from Peters (1964),
default: None, type:
Function, optional: True}Outputs:
h_0, type:
Float, description:
Strain amplitude}Implementation:
from GravitationalWaves import strain, source
import numpy as np
import astropy.units as u
# create a random collection of sources
n_values = 1500
mass1 = np.random.uniform(0, 10, n_values) * u.Msun
mass2 = np.random.uniform(0, 10, n_values) * u.Msun
lumDist = np.random.normal(8, 1.5, n_values) * u.kpc
freqOrb = 10**(-5 * np.random.power(3, n_values)) * u.Hz
ecc = 1 - np.random.power(3, n_values)
sources = source.Source(m_1=m_1, m_2=m_2, ecc=ecc, dist=dist, f_orb=f_orb, gw_lum_tol=0.05, stat_tol=1e-2,
interpolate_g=True, interpolate_sc=True})
h_c_n4 = sources.get_h_c_n(harmonics=[1, 2, 3, 4])Name: psd
Description: The psd module evaluates the effective noise power spectral density of a detector at different frequencies. The module contains the LISA sensitivity curve that can be tweaked by adjusting parameters such as the observation time, response function, and the arm length. For the Galactic confusion noise, the module uses the Robert et al. 2019 model.
Annotations:
"""
Calculates the effective power spectral density for all instruments.
"""Inputs:
f,
description:
Frequencies to evaluate sensitivity curve,
type: Float, optional:
False}instrument,
description: Instrument to use,
default: 'LISA',
type:String, optional:
True}custom_psd,
description:
Custom function for computing the PSD,
default: None, type:
Funcion, optional: True}t_obs,
description: Observation time,
default: 4, type:
Float, optional: True}L,
description: LISA arm length in metres,
default: 2.5, type:
Float, optional: True}approximate_R,
description:
Whether to approximate response function,
default: ‘False’, type:
Bool, optional: True}confusion_noise,
description: Galactic confusion noise,
default: 'robson19',
type: Various, optional:
True}Outputs:
psd, type:
Float, description:
Effective power strain spectral density}Implementation:
import GravitationalWaves.psd as psd
import numpy as np
import astropy.units as u
frequency_range = np.logspace(-5, 0, 1000) * u.Hz
PSD = psd.power_spectral_density(f = frequency_range) # get noise amplitude
noise_amplitude = np.sqrt(frequency_range * PSD)Name: utils
Description: The utils module is a collection of miscellaneous utility functions for computing conversions between variables, the chirp mass of binaries, and expressions from the Peters equation (Peters and Mathews 1963), such as the relative power of gravitational radiation and enhancement factor.
Annotations:
"""
Computes chirp mass of binaries
Compute g(n, e) from Peters and Mathews (1963) Eq.20
Converts between orbital frequency and semi-major axis using Kepler's third law
"""Inputs:
m_1,
description: Primary mass,
type: Float, optional:
False}m_2,
description: Secondary mass,
type: Float, optional:
False}n,
description: Harmonic(s) of interest,
type: Int, optional:
False}e,
description: Eccentricity,
type: Float, optional:
False}f_orb,
description: Orbital frequency,
type: Float, optional:
False}Outputs:
m_c, type:
Float, description:
Chirp mass}g, type:
Array, description:
Power of gravitational radiation at nth harmonic}f_orb, type:
Float, description:
Orbital frequency}a, type:
Float, description:
Semi-major axis}Implementation:
from GravitationalWaves import utils
import astropy.units as u
mass1 = 11 * u.Msun
mass2 = 11 * u.Msun
m_c = utils.chirp_mass(m_1 = mass1, m_2 = mass2)Name: visualizations
Description: The visualization module uses matplotlib and seaborn to plot quick visuals for analyzing of a collection of sources. The plots include 1- and 2-dimensional distributions and the LISA sensitivity curve.
Annotations:
"""
Plot a 1D distribution of `x` and `y`
Plot a 2D distribution of `x` and `y`
Plot the LISA sensitivity curve
"""Inputs:
x,
description:
Variable to plot on the x axis, type:
Float, optional: False}y,
description:
Variable to plot on the y axis, default:
None, type: Float,
optional: True}weights,
description:
Weights for each variable pair, type:
Float, optional: False}disttype,
description:
Type of distribution plot to use, type:
{ "scatter", "kde" }, optional:
False}frequency_range,
description:
Frequency values to plot sensitivity curve,
type: Float, optional:
False}y_quantity,
description: Quantity to plot on y axis,
default: None, type:
{ "ASD", "h_c" }, optional:
True}Outputs:
fig, type:
matplotlib Figure, description:
Figure of distribution}ax, type:
matplotlib Axis, description:
Axis of distribution plot}Implementation:
import GravitationalWaves.visualization as vis
import numpy as np
import astropy.units as u
fig, ax = sources.plot_source_variables(xstr = "f_orb", ystr = "snr", disttype = "kde", log_scale = (True, True), show=False, which_sources = sources.snr > 0)
right_ax = ax.twinx()
frequency_range = np.logspace(np.log10(2e-6), np.log10(2e-1), 1000) * u.Hz
vis.plot_sensitivity_curve(frequency_range = frequency_range, fig = fig, ax = right_ax)Name: wavesim
Description: The gwsim module contains functions for simulating gravitational waves from binary black holes. The procedure includes parts for simulating the inspiral portions, pre-matching and post-matching parts for simulating the merger/ringdown portions, and parts for matching together the inspiral and merger/ringdown portions.
Annotations:
"""
Parts of the procedure for simulating the inspiral portions of gravitational waves, collected into modular functions.
"""Inputs:
logMc,
description: log (Mc/Msun),
type: Float, optional:
False}q,
description: Eccentricity,
type: Float, optional:
False}flow,
description: Frequency flow,
default: 10.0, type:
Float, optional: True}merger_type,
description:
Binary black hole or neutron star merger,
type: String, optional:
False}D,
description: Distance,
default: 100.0, type:
Float, optional: True}Outputs:
i_m_amp, type:
Float, description:
Strain envelope amplitude}i_m_freq,
type: Float, description:
Frequency over time for whole waveform}With this package, users can
The most basic use case for GravitationalWaves is to calculate the
signal-to-noise ratio of a single stellar-mass binary system. We do this
by using the package’s source module to generate a Source class. First
we must import the package module. As soon as we use import statements,
we can use modules, which can be built-in, installed third-party
modules, or internal project modules. The best import statement is
simply import modu. Using
from modu import func is a way to pinpoint the function to
import and put it in the local namespace.
To import the file source.py in a directory, we use the
statement import GravitationalWaves.source. This statement
looks in the directory for __init__.py and then for
source.py, executing all top-level statements in both
files. Consequently, any variable, function, or class defined in
source.py is available in the local namespace.
import GravitationalWaves
import GravitationalWaves.source
import astropy.units as u
source1 = source.Source(m_1 = 11 * u.Msun,
m_2 = 11 * u.Msun,
ecc = 0.3,
f_orb = 1e-4 * u.Hz,
dist = 9 * u.kpc,
interpolate_g = False)
source1.get_snr()For this example, GravitationalWaves checks whether the source is eccentric/circular and evolving/stationary, and chooses the fastest method to accurately calculate the SNR.
In the next example, we use GravitationalWaves to calculate the detectability of a collection of sources. We first import the package, then create a random set of possible LISA sources.
import legwork.source as source
import numpy as np
import astropy.units as u
# create a random collection of sources
nVals = 1800
mass1 = np.random.uniform(0, 12, n_values) * u.Msun
mass2 = np.random.uniform(0, 12, n_values) * u.Msun
eccInit = 1 - np.random.power(6, n_values)
lumDist = np.random.normal(9, 1.5, n_values) * u.kpc
orbFreq = 12**(-6 * np.random.power(3, n_values)) * u.HzUsing these random sources, we can instantiate a Source class.
sources = source.Source(m_1 = mass1, m_2 = mass2, ecc = eccInit, dist = lumDist, f_orb = orbFreq)Now, we can calculate the SNR (signal-to-noise ratio) for these sources. This function splits the sources based on whether they are stationary/evolving and circular/eccentric.
snr = sources.get_snr(verbose = True)The SNR values are now stored in sources.snr, and we can
mask any values that don’t meet some detectable threshold. In the
following, we set the threshold to 7.
detectableSources = sources.snr > 7
print("{} out of the {} sources are detectable".format(
len(sources.snr[detectableSources]), nVals))Overall, this project taught me how to structure a Python project. Structuring GravitationalWaves means making decisions about how the project can best achieve its goals: simulating, analyzing, and detecting gravitational waves. In other words, the structure requires writing clean code, with clear logic and dependencies, and organizing files and folders in the file system. Writing clean code with clear logic is part of the architectural task of building out the different components of a project and their interactions. Organizing a project requires making decisions, such as which features should go into a module, how data should flow through the project, and which features and functionality should be grouped together.
By structuring the Python project, I learned how to divide the code base into clean, efficient modules using packages. This step requires building modules as an abstraction layer to divide the code into sections containing related data and functionality. For example, this project separates the layers that deal with plotting results from the layers that deal with miscellaneous utility functions by grouping all plotting functions in the visualization module and all utility functions in the utils module. Overall, the straightforward module importing model and packaging system provided in Python makes it notably easier to structure a Python package.
legwork and
riroriro package so that users can both simulate
theoretical gravitational wave forms and analyze existing sources.
Overall, GravitationalWaves is designed to help calculate binary sources of gravitational waves, whether simulated or LISA-like. This package aims to provide a way to study and better understand the detectability of such compact object binaries. Future work includes adding more functions, equations, and modules to the package to implement and analyze gravitational-wave emission, gravitational wave strain, SNR, and visualization modules, and see these implementations’ effect on orbital evolution. This will require adding more modules and functions to compute higher-level linear algebra-based equations and mathematical models.